Controllable CO2 electrocatalytic reduction via ferroelectric switching on single atom anchored In2Se3 monolayer

Efficient and selective CO2 electroreduction into chemical fuels promises to alleviate environmental pollution and energy crisis, but it relies on catalysts with controllable product selectivity and reaction path. Here, by means of first-principles calculations, we identify six ferroelectric catalysts comprising transition-metal atoms anchored on In2Se3 monolayer, whose catalytic performance can be controlled by ferroelectric switching based on adjusted d-band center and occupation of supported metal atoms. The polarization dependent activation allows effective control of the limiting potential of CO2 reduction on TM@In2Se3 (TM = Ni, Pd, Rh, Nb, and Re) as well as the reaction paths and final products on Nb@In2Se3 and Re@In2Se3. Interestingly, the ferroelectric switching can even reactivate the stuck catalytic CO2 reduction on Zr@In2Se3. The fairly low limiting potential and the unique ferroelectric controllable CO2 catalytic performance on atomically dispersed transition-metals on In2Se3 clearly distinguish them from traditional single atom catalysts, and open an avenue toward improving catalytic activity and selectivity for efficient and controllable electrochemical CO2 reduction reaction.


The intrinsic stability of α-In2Se3 and structural polymorphs
The In-Se compounds have many possible structural polymorphs, such as the phases of α, β, γ, δ, κ-In2Se3, which have been determined by X-ray diffraction and TEM. 1 However, α-In2Se3 is the ground-state phase and is stable at the room temperature from both theoretical and experimental perspectives, as elaborated below.
Theory: To examine the stabilities of different structural polymorphs, we calculated the total energies of six possible phases of the In2Se3 monolayer (Figure S1), including the β', β, α and α', zincblende, and wurtzite phases. The α-In2Se3 monolayers with ferroelectric polarizations have the lowest total energies, consistent with the recent theoretical work by Ding et al. 2 , indicating α-In2Se3 to be the most stable phase. Experiments: Different phases In2Se3 (α, β, γ, δ, κ-phase) have been experimentally synthesized, 1 but under distinct fabrication conditions. Past reports have explicitly indicated that α-In2Se3 is the room-temperature phase, 3,4 while β, γ, δ -phases are high-temperature phases. 5 Phase transformation can be achieved via the path of , 1 which also indicates that α phase is the stable room-temperature structure. This point has been further verified by the recent synthesis of α-In2Se3 layers that have been taken to fabricate ferroelectric devices. 6,7,8 Cui et al. have pointed out that the cooling rate is critical S3 for obtaining the α phase, which is stable at room temperature. 9 Therefore, α-In2Se3 is rather stable at room temperature and can be well prepared by controlling the synthesizing temperature.

Total energy of possible In/Se phases in the presence of metal adatoms and in the electrochemical environment
In order to check if α-In2Se3 is still the global minimum energy phase in the presence of metal adatoms and in the electrochemical environment, we calculated the total energies of other four phases (β', β, wurtzite, and zincblende) with presences of Pd, H* or OH*. As shown in Figure S2, for all the studied cases, the α-In2Se3 has the lowest total energy among these five phases. For other phases with different In/Se ratio, we only considered the InSe monolayer since it was experimentally synthesized. 10,11 Due to the different species and number of atoms, we cannot simply use the total energies to make a comparison between Pd@InSe and Pd@α-In2Se3. Herein the cohesive energy is employed instead.
All these results above indicate that, α-In2Se3 is the global minimum energy phase, even in the presence of adatoms and in electrochemical environment.

Stability of α-In2Se3 under harsh environments
We have evaluated the electrochemical stabilities of α-In2Se3 monolayer by the dissolution potential Udiss, 12,13,14 which is defined as where 9:44 。 ( ) and n are the standard dissolution potential of In/Se bulk and the number of electrons involved in the dissolution, respectively, which can be obtained from the NIST database. 15 ! is the formation energy of In/Se atom in the In2Se3 monolayer, given by: where 23+<=>? , 45+<=>? are the respective total energies of the In and Se atoms in their most stable bulk structures, 23 ! ;5 " is the total energy of the In2Se3 monolayer. According to this definition, materials with 9:44 > 0V vs SHE are regarded as electrochemically stable under acidic conditions. The 9:44 of both In and Se in the In2Se3 monolayer are positive (see Table   S3), indicating the electrochemical stability of the In2Se3 monolayer.  Figure S3.

Figure S3
Comparison of CO2RR pathways on α and α' phases of Pd@In2Se3.

Figure S5
These 29 kinds of transition metals were chosen to modify the In2Se3 surface to increase its reactivity.

Metal substituted In2Se3
To determine possible metal substation on the In sites, we have simulated the formation of In vacancies. The formation energy of In vacancy and the diffusion barriers of an In atom removal from In2Se3 (see Figure S6)   In addition, we find that the metal substituted In2Se3 is not suitable for the electrocatalytic CO2RR, even if we do not consider the difficulties of In vacancy formation and metal substitutions. According to the free-energy profile for CO2 electrochemical reduction reactions along the minimum energy path at 0 V (vs. RHE) on Pd substituted P↑-and P↓-In2Se3 (namely Pd substitutes the In vacancy), the energy barriers of CO2 + *→ OCOH* are S10 up to 1.17 and 1.88 eV on Pd-doped P↑-and P↓-In2Se3, respectively, due to the full occupations of Pd-d orbitals (see Figure S7), leading to the low CO2RR activities of these catalysts. Overall, we conclude from our calculations and analysis that (i) In vacancies in In2Se3 are hard to form due to the high formation energies and diffusion barriers, making metal substituted In2Se3 rare in reality and (ii) metal substituted In2Se3 is unsuitable for CO2RR.
Therefore, we focus our discussions in this work on metal anchored In2Se3 on the surface.

Figure S8
Side and top and views of optimized configuration of TM@P↓-In2Se3 (a, c, e, g, and i) and TM@P↑-In2Se3 (b, d, f, h, and j). The 3D differential charge density plots of TM@ In2Se3 are also shown in this figure, which are obtained with these optimized configurations.
OPQR * is defined as: where 6PQR6 and T ! are total energy of the hexafluoroacetylacetonate and hydrogen molecules.
The calculated values of ∆< are -0.70 and -0.15 eV for Pd@P↓-In2Se3 and Pd@P↑-In2Se3, respectively. We can thus conclude that: 1) the adsorption of Pd atom on In2Se3 is more stable, than that of Pd(hfac)2; 2), the Pd@In2Se3 could be synthesized by using Pd(hfac)2 as the molecular precursors.
Using the same methods, we also calculated the energy differences ∆< for Rh and Ni based catalysts, see Table S4, all of them are negative, indicating the stability and possible synthesis based on the molecular precursors.
For Re, Zr and Nb, the corresponding molecular precursors for single atom catalysts synthesis are not found from the previous literatures, but we noticed that the metal rods were used as the sources for the fabrications of single atom catalysts, 17 we therefore keep the binding energies of these SACs calculated from isolated atoms in our work, which can be used to reflect the binding strengths. S14

The clustering energy
The clustering tendency of Pd atoms on the surface is estimated by the clustering 20,21,22 which is defined as the difference in binding energies between a single metal atom ( <,4:3 ) and the metal dimer ( <, 9:] ): where ( <,4:3 ) and ( <, 9:] ) are defined as: where ^9+<=>? is the chemical potential of the Pd atoms in their bulk phase, and ^9 ! /23 ! ;5 " represents the total energy of the substrate with a Pd dimer. According to the definitions, negative values of Z>=4[5\ mean that the metal cluster does not tend to form. The calculated values of Z>=4[5\ are -0.3 eV and -0.04 eV for Pd@P↓-In2Se3 and Pd@P↑-In2Se3, respectively. To further analyze the stability of single-atom adsorption, we have performed first-principles finite-temperature molecular dynamics simulations of two dispersedly adsorbed Pd atoms on P↓-In2Se3 (2Pd@P↓-In2Se3) with a Nose-Hoover thermostat at 300 K.
The fluctuations of the temperature and the total energy as a function of the simulation time are given in Figure S11b. The two Pd atoms maintain dispersedly adsorbed features for at least 15 ps. The distance between the two Pd atoms stays essentially unchanged. All these results confirm the dynamic stability of the single-atom adsorption at room temperature.

S17
The kinetic Monte Carlo (kMC) simulations KMC simulation method. The details of kMC method have been described in previous works. 23,24,25 In this case, we start with a perfect α-In2Se3 monolayer with evenly distributed single metal atoms. The simulation size of the α-In2Se3 monolayer is ~400×400 Å 2 with the periodic boundary conditions and the density of Pd monomers is ~4% (the same model and method also used for Nb case), as shown in will not agglomerate at room temperature (i. e. 293K) environments for at least 100 seconds.

Table S5
The values of activation barriers used in the KMC simulations, which are obtained from the DFT calculations: binding energies of metal monomers and metal atoms in metal clusters ( : in eV/atom) and diffusion barriers of metal monomers and metal atoms in metal clusters ( : in eV/atom).

Pd@P↓-In2Se3
Pd monomer - The kMC result for Pd@P↑-In2Se3 is similar (Figure S13), although the starting temperature to form the Pd cluster is much lower than that of Pd@P↑-In2Se3 due to lower diffusion barrier (see Table S5), no Pd clusters form at T≤300K for 100 seconds, indicating the stability at the room temperature environments. Since Pd@P↑-In2Se3 has the lowest migration barrier of 0.84 eV (see Table 1

Local metallization in TM@In2Se3
Since ferroelectricity generally stems from the offset positive/negative charge centers, the switchability of polarization is closely related with the electronic structures, and the semiconducting or insulating feature is essential for the ferroelectrics.
The catalysts studied here are locally metallized, with the metallic states mainly from the anchored metal atoms and their immediate surrounding area. To prove this point, we have built a large supercell with Pd anchored 6×6 P↓-In2Se3 monolayer as shown in Figure S17 (left). Based on the projected DOS analysis (see the right of Figure S17

Local vs. whole ferroelectric switching
To check if the local or whole ferroelectric switching is energetically preferred, we have built a supercell with HOCHO* adsorbed Zr@6×6 P↑-In2Se3 ( Figure S22a) and calculated the energy differences of the structure with ferroelectric switching at only small ( Figure S22b) or large (Figure S22c) local area around Zr anchored site and (Figure S22d) entire lattice.
The calculated results indicate that the local ferroelectric switching is indeed energetically prohibitive, and the cases (b) & (c) have much higher energies than that of (d).
However, the ferroelectric switching of the entire lattice is not difficult, even when the single metal atoms are fairly apart from each other (24.41 Å): the energy difference is 0.4 eV/6×6cell or 0.045 eV/2×2cell, it is comparable to the data shown in the manuscript.
Therefore, we conclude that the ferroelectric switching can readily occur throughout the entire lattice due to the energetic preference when HOCHO* exists.
where 4=< and ef+4=< are the total energies of the substrate without and with the added TM single atom, respectively. ef is the total energy of the isolated TM atom. The binding energy for the adsorbed CO2 molecule (Eb-CO2) onto TM@In2Se3 (TM=Nb, Ni, Pd, Rh, Re and Zr) is evaluated in the similar way.

Charge density difference
The charge density difference for TM@In2Se3 (TM=Nb, Ni, Pd, Rh, Re and Zr) (∆ * ) is calculated by the following equation: where 4=< and ef+4=< are the charge density of the substrate without and with the added TM single atom, respectively. ef is the charge density of the added TM atom. The charge density difference for the CO2 molecule adsorbed TM@In2Se3 (TM=Nb, Ni, Pd, Rh, Re and Zr) (∆ 0 ) is evaluated in the similar way.

d band centre
In order to get the individual orbital components, we employed the code named "splitdos.dos" to process the output file of density of states (DOS). The average d-band (d band center) shifts are calculated for the surface metal atoms for both the total d partial DOS and the orbital-resolved d partial DOS. The d band center ( 9 ) is calculated as: 39 Where is the electronic energy of states, and 9 ( ) is the electronic density of states.

Reaction free energy change
The reaction free energy change (ΔG) for each elementary step is calculated based on the computational hydrogen electrode (CHE) model proposed by Nørskov and co-workers by the following equation: 40, 41 where ∆ is the reaction energy obtained from DFT calculations. ∆ i^j and ∆ are the contributions of the zeropoint energy and entropy to ∆ , respectively, which are obtained from the vibrational frequency. T is the temperature and taken as 298.15 K, and ΔS is the entropy change.
The vibrational frequencies of molecules in the gas phase are taken from the NIST database, 15 and others are calculated with considering solvent effect. Tables S7-S19 present i^j and TS (at 298.15 K) of the free molecules and the adsorbed species along the most favourable reaction pathway for the TM@In2Se3 (TM =Nb, Ni, Pd, Re, Rh, and Zr) catalyst including TM-down and TM-up configurations. e and U are the number of electrons transferred and the electrode potential applied, respectively. ∆ kT is the free energy correction of pH, which is calculated as follows: where c is the Boltzmann constant, and pH = 0 is assumed in an acidic medium in this study.
The value of limiting potential ( > ) is determined by the potential-determining step (PDS), which has the most positive ∆ (∆ ]lA ), as computed as follows:

Elementary steps in CO2RR and HER
All the hydrogenation reactions (R1~R29) considered in the search process for the minimum energy reaction pathways of the CO2 reduction reactions (see Figure S21) can be written as: The hydrogenation reactions (R30~R31) for HER can be written as: Therefore, when U = 0 V and pH = 0, ΔG for each elementary step (R1-R31) can be rewritten as: The Gibbs free energy (G) of the adsorbates on the TM@In2Se3 catalysts along the minimum energy reaction pathway for both CO2 RR and HER are listed in Table S8~S19. No 1 --"Yes" means these metal atoms match the criteria. "No" means these metal atoms mismatch the criteria. Specifically, "No 1 " means these metal atoms desorb from the In2Se3 surface at both up and down polarization phases; "No 2 " means these metal atoms could stably adsorb on the In2Se3 surface at only one polarization phase; "No 3 " means the favourable adsorption sites of these metal atoms obviously change after the polarization switch; "No 4 " means these metal atoms can not activate CO2 molecules at both up and down polarization phases; "No 5 " means these metal atoms can activate CO2 molecules at only one polarization phases, but the species, composed by the activated CO2 molecules and the metal atoms, desorb from the In2Se3 surface. "--" means these metal atoms are not considered for the second screening procedure, because they mismatch Criteria 1.

The number of catalytic active sites
The simulations were conducted to study the maximum active sites of Pd adatoms on a 4×4 In2Se3 supercell. Since the centers of the six-membered ring have been identified as the energetically preferred dopant sites (Figure 1), where the metal adsorbents are uniformly dispersed (rather than forming cluster) as the catalytic active sites. Based on the results, the formation energies are calculated to figure out the maximum active sites of SAC, with Pd adatoms evenly distributed 4×4 In2Se3 supercell as the representative example. As shown in Additional Pd atom on the surface will raise the formation energy due to the reduced distance between the nearest neighbors and the Coulomb repulsion. For example, the fifth Pd adsorbent on 4×4 In2Se3 leads to the formation energies of -1.65 and -1.12 eV/Pd atom for two polarization states. We therefore conclude that In2Se3 surface can host a large number of catalytic active sites, which are uniformly dispersed with the converge up to 25%.

The dependence of catalytic activity on metal atom concentrations
To investigate the dependence of catalytic activity on metal atom concentrations, we have calculated the free-energy profile for CO2 electrochemical reduction reactions along the minimum energy path at 0 V (vs. RHE) on Pd@In2Se3 by using 2×2, 4×4, and 6×6 α-In2Se3 S44 supercells. As shown in Figure S24, the difference of the free-energy profiles for the CO2RR between Pd@ In2Se3 (4×4) and Pd@In2Se3 (6×6) are within 0.05 eV, and the converged freeenergy profiles indicate that the catalytic activity is preserved when the doping concentration is in the dilute limit. On the other hand, while the binding strengths of the reaction intermediates for the dilute doping concentration are stronger than those for the high concentration doping case, the rate limited steps for CO2RR are the same, and the overpotential difference of CO2RR between the high and low concentrations is less than 0.13 eV. These results confirm that the CO2RR activities of metal doped In2Se3 with high and low concentrations are quite similar, and the doping concentration on metal adsorbed In2Se3 has little effect on our main conclusions.

PBE vs. PBE+U
The inclusion of the Hubbard-U term via, e.g., the DFT+U approach, may be more suitable for systems with highly localized orbitals. However, the DFT+U approach also suffers from a strong (linear) dependence of the energetics on the choice of the value of the parameter U, and on the choice of the localized projector functions that enter the definition of the U-dependent energy term. For example, the reduction energy (ΔH) of CeO2→Ce2O3 process can vary between −5.1 (U = 0 eV) and −1.9 eV (U = 5.0 eV) using the DFT+U method 33 , while the GGA-PBE value of −4.18 eV is in good agreement with the experimental measurements (−3.57 to −4.03 eV). On the other hand, the U-value is usually chosen based on its accuracy in reproducing the electronic structures (i.e., experimental band gap) of the bulk materials. However, to simulate catalysts, it is better to choose U to fit the energy of the oxidation-reduction, since catalytic processes are controlled by energy difference 34 . The specific case in this work, namely CO2RR on SAC surfaces, involves complex surfaceadsorbate interactions, under which bulk-property derived U values in a locally changing surface environment may not adequately describe the reaction energetics 35,36 .
Note that the results based on the GGA-PBE (the method used in this work) showed very good performance in understanding the reaction mechanisms and activity trends observed in experiments 37,38 . To verify the accuracy of the PBE results, we also investigated the CO2RR pathways on the Pd@In2Se3 with the PBE+U method in which the previously validated U value of 8.00 eV was employed for the Pd 4d orbital 42 . Although the PBE+U results of limiting potentials are slightly (less than 0.2 V) larger than the PBE ones, the computed theoretical final product, reaction path, potential-liming step as well as the variation of limiting potential caused by polarization conversion are the same (see Figure   S25). Thus, the standard PBE calculations provided accurate predictions to describe the CO2RR activity.

Figure S25
Comparison of the CO2RR pathways on Pd@In2Se3 by using the PBE and PBE+U methods.

PBE vs. RPBE
It has been reported that the RPBE functional suggested by B. Hammer et al. could improve the chemisorption energetics of atoms and molecules on transition-metal surfaces. 43 In the present work, to evaluate the accuracy of the PBE functional, we also investigated the CO2RR pathways on the Pd@In2Se3 with the RPBE functional. As shown in the Figure S26, except for the slightly increased limiting potential (less than 0.11 V), RPBE gives the same results as PBE, including the final product, reaction path, potential-liming step, as well as the variation of limiting potential caused by polarization conversion. Therefore, the conclusions obtained with the PBE functional hold true when the more accurate RPBE functional is used.

Figure S26
Comparison of CO2RR pathways on Pd@In2Se3 by using the PBE and RPBE methods.